Analytical approach to directed sandpile models on the Apollonian network 
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We investigate a set of directed sandpile models on the Apollonian network, which are inspired on the work 
by Dhar and Ramaswamy (PRL 63, 1659 (1989)) for Euclidian lattices. They are characterized by a single 
parameter q, that restricts the number of neighbors receiving grains from a toppling node. Due to the geometry 
of the network, two and three point correlation functions are amenable to exact treatment, leading to analytical 
results for the avalanche distributions in the limit of an infinite system, for q = 1,2. The exact recurrence 
expressions for the correlation functions are numerically iterated to obtain results for finite size systems, when 
larger values of q are considered. Finally, a detailed description of the local flux properties is provided by a 
multifractal scaling analysis. 



I. INTRODUCTION 

The interaction networks of many real systems with large 
number of basic units are often found to display power-law 
distribution of node degrees and small- world property fHO]. 
Examples stem from most different areas, as electric power 
distribution, food webs in ecology, information flow in the 
internet, interaction among financial institutions, and so on 
JH0]. In recent years, complex networks have also attracted 
attention as alternative topological structures to ordered eu- 
clidian lattices, on which many physical models can be de- 
fined. These structures offer a suitable scenario to mimic the 
effect of geometry in real systems, and have already been used 
in the investigation of the properties of magnetic JHIaH] and 
electron [8] systems. 

Understanding the stability of complex networks becomes 
of relevance for the management of natural and human built 
systems, as it can provide guidelines to avoid an irreversible 
collapse and to enhance the robustness of their structure. An- 
other issue that deserves attention is the occurrence of events 
that may cause permanent or temporary damages on the net- 
work, which can be interpreted as avalanches within the pro- 
posed self-organized criticality (SOC) scenario |@]. It is well 
known that a typical signature of SOC systems is the possi- 
bility of occurrence of a very large avalanche that can extend 
itself over the whole network, causing its breakdown. Spe- 
cific sandpile models defined on complex networks have been 
recently investigated [10], as well as models where the net- 
work is not fixed, but the set of connections evolves slowly 
with time Hill . In the first case, avalanches refer to the mo- 
tion of mass units from one node to its neighbors, while in the 
last approach, avalanches refer to bursts of rewiring connec- 
tions among the network nodes. It is noteworthy the recent 
attempts to use SOC concepts with respect to brain activity, 
both in Euclidian and scale-free networks Il2lll3lll4ll . 

It is well known that direct models, like the one proposed by 
Dhar and Ramaswamy lfl5ll . constitute one of the few classes 
of SOC models that can be exactly solved on Euclidian lat- 
tices. This is essentially related to their Abelian property, ac- 
cording to which the effect of two successive grain additions 



on the lattice does not depend on the order. In the context 
of complex networks, the Apollonian packing problem lfl6ll 
inspired the introduction of the so-called Apollonian network 
B17lll8ll . Besides displaying both scale free and small-world 
features, the hierarchical geometry of this network enables the 
derivation of tractable analytical expressions for a variety of 
equilibrium and dynamical models [19]. This leads either to 
exact results or to recurrence relations that can be numerically 
iterated. 

In this work, we analyze the avalanches of directed sandpile 
models on the Apollonian network. We make use of properties 
of these specific network and model to derive, in a first place, a 
series of exact results for the distribution of avalanches. Then, 
these results can be extended, with the help of the numeri- 
cal iteration of the obtained recurrence relations, to illustrate 
more general situations. More precisely, we are able to in- 
vestigate the fine details of the local mass flux, deriving the 
appropriate multifractal spectra that describe the scaling prop- 
erties of the flux. 

This work is organized as follows: In Section II we intro- 
duce our model, discussing the role played by the number of 
levels q, in the Apollonian hierarchy, that limit which nodes 
can receive mass from a toppling neighbor. We also derive 
the basic expressions for the two and three-point correlation 
functions that allow for the derivation of the local and total 
fluxes. Results for the total flux, obtained by numerical it- 
eration, are discussed in Section III, for 1 < q < 6 . They 
are then compared with analytical expressions derived for the 
q = 1 and q = 2. In Section IV, a multifractal approach is 
used to present the scaling properties of the flux for the dis- 
tinct values of q. Finally, Section V closes the paper with our 
concluding remarks. 



II. DIRECTED SANDPILE MODELS ON THE 
APOLLONIAN NETWORK 

The planar Apollonian network 11711 is obtained from the 
classical Apollonian space-filling packing of circles l2(ill . by 
associating nodes with the centers of the circles, and draw- 
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ing edges between nodes corresponding to pairs of touching 
circles. This iterative building process is illustrated in Fig. Q] 

The directed sandpile model of Dhar and Ramaswamy 01511 
associates with each site x of a hypercubic lattice a height 
variable z(x), which is increased by 1 when a grain is added 
to x. If z(x) exceeds a critical value Zc, the site topples, and 
the height variables at the £ nearest neighbors of x along a 
preferred direction increase by 1, while z(x) decreases by I. 
Without loss of generality, z c can be chosen to equal I, The 
existence of a preferred direction is essential to the exact solv- 
ability of the model, not only in its original form but also in 
generalized versions 12 ill . 

In the Apollonian network, the building process offers an 
obvious choice of a preferred direction. We define the nth 
layer of the network as the set of sites added in the nth it- 
eration of the process, and we postulate that, when a site at 
a given layer topples, only sites in subsequent layers can re- 
ceive grains. However, the Apollonian network has the pecu- 
liar property that each site in a given layer is connected to at 
least one site in each subsequent layer. Thus, in the thermo- 
dynamic limit, any site has an infinite number of neighbors 
in subsequent layers, leading to an infinite critical height. In 
order to obtain a finite value of z c , we impose the restriction 
that only neighbors in the first q subsequent layers can receive 
grains when a site topples, the remaining connections being 
inactive; see Fig. [2] This leads to a g-dependent value of the 
critical height z c , which is the same for all sites in the net- 
work, provided we forbid the addition of grains to the sites in 
the original triangle (layer n = 0). (For q from 1 to 6, we have 
Z c =3, 9, 21, 45, 93, and 189.) Thus, all allowed sites have an 
equivalent set of neighbors in their subsequent layers, and we 
can study the properties of avalanches by choosing any ref- 
erence site xo. For convenience, we choose xo to be the site 
located at the geometrical center of the network (layer n = 1). 

As in the original directed sandpile model, we define a 
two-point correlation function Go (x;xo) which measures the 
probability that a site x topples in the SOC state due to an 
avalanche originated by adding a grain at xo. Since the prob- 
ability that a site topples, provided that r of its backwards 
neighbors have toppled, is equal to r/z c , Go obeys the recur- 
sion equation 



Go (x;x ) 



1 

Zc 



Y,' G o (y;xo) + 5 XjX 



(i) 



with the primed summation running over all sites from which 
x can receive grains, according to the g-layer rule. Since 



G (x ;x ) 



1 

Zc' 



(2) 



the existence of a preferred direction allows us to solve Eq. 
((TJ for all Go (x;xo), at least numerically. 
The flux through the nth layer is given by 



<K«) = £g (x;x ). 



(3) 



Contrary to what is observed in hypercubic lattices, here (n) 
generally depends on n, although it becomes asymptotically 



constant for n> 1, as we show below by numerical and ana- 
lytical calculations. If m (n) is the average number of sites in 
the nth layer that topple when at least one of them does, we 
can write 



<|)(n) = m (n) p (n) . 



(4) 



in which p (n) is the probability that, in the SOC state, an 
avalanche started at xo reaches layer n. 
If we assume that 



p{n)~n a , 



(5) 



with some exponent a, the asymptotic constancy of (j) (n) al- 
lows us to conclude that 



m (n) 



1 



p(n) 



(6) 



Thus, the average mass of an avalanche reaching n layers 
scales as 



M (n) = £m (?) ~ f dt t a - n a+1 , 



(7) 



and the probability that the total mass of an avalanche exceeds 
M can be written as 



p(M) =p(n(M))~M~ 



(8) 



Finally, we obtain for p(M), the probability distribution of 
avalanches with size M, 



p(M) 



dp(M) 
dM 



■M ¥Fo =M x 



(9) 



The exponent a can be calculated from the mean square 
flux 



4>(n) = [m(n)f p(n) — n a , 



(10) 



which is related to the three-point correlation function 
G(xi,X2;xo), defined as the probability that sites Xi and X2, 
both in the same layer, topple due to an avalanche started by 
adding a grain at site Xq. Explicitly, we have 



4>(n)= £ G(xi,x 2 ;x ). 

xi ,X26n 



(ID 



As in the case of Go (x;xo), we can write for G(xi,X2;xo) 
a recursion equation, 

G(xi,x 2 ;x ) = ^ ^ 'G(yi,y 2 ;xo), (12) 

with the primed summation running over all sites from which 
xi or X2 can receive grains, according to t he g-lay er rule. This 
last equation can be solved by the Ansatz 11221,12311 

G(xi,x 2 ;x ) =£/(y;xo)Go(xi;y)Go(x 2 ;y), (13) 
y 



Figure 1 : Building process of the Apollonian network. 



with the function / (y;xo) determined by the condition 

G(x,x;x ) =Go(x;x ), (14) 

which leads to 

^/(y;xo)G (x;y)G (x;y) = G (x;x ). (15) 
y 

Summing over all sites x in the same layer n, using Eq. (01 
and the fact that G (x; y) = G (x — y + Xrj; Xo), we can rewrite 
Eq. (|T5j as 

£f(0#(«-;+i) = <K«), (16) 

1=1 

in which 

^W=£/(y) and ^(f) = ^Go(x;xo)Go(x;x ). 

yer xet 

(17) 

Starting from n = 1, Eq. ( [ToT l can be solved recursively for 
F (n). By substituting Eq. ( fT3l into Eq. ( fTTt . we can express 
<I> (n) in terms of F (n), 

<S>{n) = j^F{t)[$(n-t+\)] 2 . (18) 

r=l 

The scaling behavior of <I> (n) determines the exponent a. 

The case q = 1 is immediately solved. In this limit, the 
Apollonian network (with the three original vertices removed) 
reduces to a Cayley tree with coordination number equal to 4, 
as shown in Fig. [2] The two-point correlation is easily seen to 
satisfy 

Go(x;x ) = — , xen, (19) 
so that the average flux is (j)(n) = 1 /3, Vn, leading to 

*(") = ^TT< F(l)=3,F(n)=2(n>l). (20) 
Thus, the mean square flux is given by 

«>(") = ~ + ~", (2D 

corresponding to a = 1 (x = 3/2), characteristic of the mean- 
field behavior associated with the directed model in Bravais 
lattices with dimension d>4. For the purpose of comparison, 
the corresponding exact values in d — 2 are a = l/2,T = 4/3. 

In the next two sections, we discuss the properties of the 
model for q > 1 . 



it — 0, type 1 

O " = 1, type 1 

• n = 2, type 1 

• n = 3, type 1 

• ii = 3, type 2 




Figure 2: Apollonian network with q = 2. Dotted lines correspond to 
inactive connections; thick lines indicate connections between sites 
in adjacent layers, while dashed lines connect sites separated by 2 
layers. Note that there are two types of sites in layer n = 3. If q — 1, 
also dashed lines become inactive. 



III. AVERAGE BEHAVIOR 

For q>2, sites in the same layer are no longer equivalent, 
since we are preserving the underlying topology of the Apol- 
lonian network as defined by the building rule. Instead, those 
sites are naturally grouped in different classes, defined by the 
structure of their connection to sites in previous layers. In 
principle, this makes the model amenable to analytical treat- 
ment. As we show in Appendix lAl the analysis for q = 2 is 
already somewhat intricate, but it lends support to a series of 
conclusions we obtain by numerical calculations. These are 
performed by building an Apollonian network with up to 16 
layers (corresponding to 21 523 363 sites), imposing the q- 
layer rule, and solving recursively Eqs. (Q]i and ( fT6b . From 
this, we can calculate both the mean flux (]) (n) and the mean 
square flux <t> («) as functions of the layer index n. (In Sec. [TV] 
we study the local properties of the flux.) 

The first conclusion to emerge from our numerical analysis 
is that the mean flux (j) («) becomes asymptotically constant 
for large n, as already mentioned in Sec. [TT] This is evident 
in Fig. |3(a)| where we plot, for several values of q, the ratio 
between (]) (n) and the corresponding (constant) result for q — 
1. Notice the oscillations in (j) (n) for small values of n. These 
are related to the fact that, as the number of neighbors of a 
site in a given subsequent layer increases with the layer index, 
so does the fraction of grains received by each layer when the 
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(a) 




(b) 



Figure 3: (a) Mean flux as a function of the layer index n, for different 
values of q, divided by the mean flux for q = 1 . (b) Corresponding 
curves for the mean square flux. 



site topples. Since we have a single site through which grains 
enter the system (in layer n — 1), we note that, for q > 1, the 
mean flux reaching layer n = 2 drops in comparison with the 
total flux, then increases from n = 2 to n = q + l, dropping 
again for n = q + 2. But this second drop is smaller, because 
sites in that layer receive grains from all q previous layers. As 
a result of this process, the oscillations are smoothed out for 
sufficiently high values of n. 

Corresponding curves for the mean square flux are shown 
in Fig. |3(b)| Again the curves oscillate for small values of 
the layer index n, but approach a constant value for large n, 
showing that <t>(n) always satisfies the scaling form 



*(n). 



(22) 



This means that asymptotically the distribution of avalanche 
sizes follows a power-law with exponent a = 1, irrespective 
of q. 

This prevalence of a mean-field behavior could be antici- 
pated on the basis of the tree-like topology of the lattice ob- 
tained by imposing the g-layer rule. A similar situation arises 



in other sandpile models on different forms of decorated Cay- 
ley trees 11231 12411 ; since the correlation length is infinite in 
the SOC state, the mean-field behavior characteristic of the 
ordinary Cayley tree (or more precisely the Bethe lattice) is 
recovered. 

In the Apollonian network, the g-layer rule defines a typi- 
cal length beyond which average properties become indistin- 
guishable from those of the model on a Cayley tree. However, 
at smaller scales, hints of the behavior corresponding to the 
genuine Apollonian network do appear, for instance in the os- 
cillations observed in the mean-flux curves. As clearly shown 
in Fig. |3(a)| («) depends exponentially on n between n = 2 
and n = q+l, 



0(n) =Ae a 



(23) 



with a ^-dependent prefactor A, but a nearly constant value of 
a ~ 0.7. The prefactor A decreases exponentially with q, since 
it is related to the inverse threshold height 1 jz c - The exponen- 
tial (rather than linear) dependence of («) is a consequence 
of the exponential increase in the number of neighbors as a 
function of the layer separation. In the q — ► °° limit, the cen- 
tral site topples only after the addition of an enormous number 
of grains, most of which are then received by sites in very dis- 
tant layers. As a consequence, all avalanches have arbitrarily 
large range. 



IV. MULTIFRACTAL PROPERTIES OF THE FLUX 

Although the average behavior of the flux reproduces that 
of the mean-field limit, the local-flux distribution reveals in- 
teresting properties already for q — 2. In Fig. |4(a)| we plot 
histograms of the local flux for n — 22. Note that, due to a 
precise identification of the distinct types of sites for the q = 2 
model, we were able to consider much larger number of nodes 
(> 10 10 ) than for the results reported in the previous section. 
The fluxes are re-scaled by the corresponding maximum flux 
0max in that layer. 

In analogy with studies on the distribution of currents in the 
incipient infinite cluster of a random-resistor network [25], it 
is interesting to evaluate the moments of the flux, in order to 
better reveal the scaling properties hidden in Fig. |4(a)| So we 
use the definition 



M k (n) 



(24) 



in which the summation runs over all sites x in the nth layer 
of the Apollonian network, 



X = G (x;x ) 



(25) 



and 0o is the initial flux. It turns out that, for all real values of 
k, the moments satisfy scaling relations given by 



M k («) 



»k" 



(26) 



with well defined coefficients u k , so that, in terms of the sys- 
tem size 



371+1 



(27) 
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Figure 4: (a) Flux distributions in a given layer, re-scaled by the 
corresponding maximum flux, for q = 2 and n = 22. (b) Moments of 
the flux as a function of n, for q = 3 and several values of k, with the 
corresponding exponents y k . Notice that there is no linear relation 
between the exponents. 




Figure 5: Plots of /(oc) for different values of q. 



that, according to Eq. d29l ), the maximum of / (6c) occurs for 
the value of a associated with k = 0, for which / (6Cfc = o) = yo- 
Indeed, for q = 1 (not shown in Fig. [5), the curve consists of a 
single point, (ECjfc = o,yo) = (1, 1), corresponding to amonofrac- 
tal behavior. Within numerical errors, that point is the maxi- 
mum of all curves, in agreement with the fact that yo = 1 for 
all values of q. 

For q>2, the left (right) end of the curves reflects the scal- 
ing behavior of the set of points associated with the largest 
(smallest) fluxes. Although not visible in the plots, the den- 
sity of points is much larger near the ends of the curves, with 
intermediate points coming mostly from values of k between 
— 1 and 0. The width of the curves increases with q, presum- 
ably diverging as q — > °°. This is related to the fact that larger 
values of q lead to a larger range of values of the flux in each 
layer of the lattice. 



V. CONCLUSIONS 



we have 

M k {L)~Ifl, (28) 

with yt — M,t/ In 3. For q = 1, all sites in a given layer n 
have the same flux 3~", so that the exponents yk are given 
by yic — 1 — k. For q > 2, on the other hand, we see from our 
numerical calculations that there is no simple linear relation 
between the exponents yk, suggesting that no single number 
characterizes the current distribution. This is a signature of 
multifractal behavior. A plot of M% (n) for q = 3 and several 
values of k is shown in Fig. |4(b)| 

To further investigate the multifractal properties of the flux 
distribution, we evaluate the multifractal spectrum /(oc), de- 
fined by a Legendre transform of the exponents y*, 

f{a)=y k +ka, H = ~- (29) 

dk 

Plots of / (6c) for several values of q are shown in Fig. [5] Note 



In this work we investigate directed sandpile models on 
the Apollonian network. Exact results were obtained for the 
avalanche distributions when q = 1 and 2. These correspond 
to the situations where an unstable node topples only to neigh- 
bors introduced in the first and the second subsequent genera- 
tions, respectively. The avalanche distributions follow power- 
law behavior, with typical mean field exponents. The iteration 
of the exact expression for the two and three point correlation 
function provides evidences for the same asymptotic behav- 
ior, regardless of the finite value of q. On the other hand, 
our results also show the emergence of large oscillatory de- 
viations due to finite size effects. This general behavior can 
be explained by noting that any finite value of q asymptot- 
ically constraints the sandpile model on the Apollonian net- 
work to the structure of a tree, where it behaves like a mean 
field model. 

The investigation of the local properties of the fluxes 
through each node shows that the network geometry induces 
a large degree of inhomogeneity in the sandpile model. This 
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effect has not been observed for the same model in Euclidian 
lattices. Nevertheless, this dependence can be accurately ac- 
counted for by a multifractal analysis. Only for q — I, when 
the model is equivalent to that defined on a Cayley tree, all 
nodes become indistinguishable, and the scaling analysis re- 
duces to a single point. 

Finally, the comparison with results found for another sand- 
pile model on scale-free networks [10] shows similarities, in 
the mean field behavior when all nodes share the same crit- 
ical height or the critical height depends locally on the node 
degree. 



Appendix A: ANALYTICAL TREATMENT FOR q = 2 

For q>2, the sites in each layer of the Apollonian network 
can be grouped in types, according to how they are connected 
to their backwards neighbors. This feature can be exploited 
in order to obtain analytical results for the behavior of the 
directed sandpile model. Here we deal with the case q = 2, 
which allows us to check our numerical results in a reason- 
ably simple way. It is clear that the treatment can be extended 
to higher values of q, with basically the same results, but a 
considerably larger amount of work. 

Layer n — 1 of the network contains only one site, while the 
three sites in layer n = 2 are all equivalent. However, already 
for n = 3 two types of sites are present: sites of type 1 receive 
grains from sites in the two previous layers, while sites of type 
2 receive grains only from the latter layer; see Fig. [2] For 
n = 4, two additional types of sites would appear, since it is 
possible that a site receives grains from sites of types 1 or 2, in 
one or two of the previous layers. It is easy to convince oneself 
that the number of site types doubles for each additional layer 
(starting at n = 2), and that the types can be labeled so that 
each site of type s has as nearest neighbors in the next layer 
two sites of type 2s — 1 and one site of type 2s. 

Denoting by g ns the value of Go (x;xo) for a site x of type 
s in layer n, and by V n>s the number of such sites, the flux 
through layer n can be written as 



,v=l 



(Al) 



In order to estimate (j) («), we must investigate the asymptotic 
behavior of both v„,. s and g ns . 

Our choice of labels allows us to write, for s = 2j — 1 (j = 1, 
2,3,...), 



and Vi.i . However, analytical results can be derived from the 
observation that g n \ behaves as 



in.\ 



(A4) 



with C = (1 + v37) / 18 - °- 393 bein 8 determined from the 
solution of the equation 



S 2 = ~(S+i). 

Consequently, g ns satisfies 



(A5) 



(A6) 



with constant prefactors A s . Moreover, the multiplicities V ns 
are such that v„,i = 3 • 2"~ 2 (n > 2) and the ratios f s = v„. 4 -/v„.i 
satisfy 



f _ v n.2j-l _ 2V„J _ f 

/2 -M - v„., - 2v,„ 
f. . - Mi = = I f . 

nj - v„, _ 2v„, _ 2JJ- 



(A7) 



We can rewrite Eq. ( IAU as 

2"- 2 b-2 

4> («) = v„,i £ f s g n , s = v„,i £ r n , m , (A8) 

s=l m=0 



with 



r„,Q = g„,i and r n , m = E f s % n > s ( m ^!)- <A9) 

.s=l+2'"-' 

Making use of the definition of F n m and of Eqs. ( 1A2| |. d A3t 
and (IA7K we can obtain the recursion equation 



r«oi — (Fn— 2.m— 2 F n — [m— 1 ) 
O 



(A10) 



Keeping in mind Eq. iA6i . we expect that T n ,„ takes the 
asymptotic form 



(All) 



Substituting this last expression into Eq. dAlOt , we con- 
clude that the constants y m satisfy the equation 



YmC - 2Ym-lC- Tlm-2 = 0, 



(A12) 



?n,2j-1 = g [ Mn 1 i 



'n-2. 



V„,2;-l = 2V„_l J, 



(A2) 

in which [w\ denotes the integer part of the number w, and, 
for s = 2 j, 



1 

?n,2j = ggn-l,j, V„, 2 > = V„_i j. 



Equations ( IA2I ) and jA3t . being recursive expressions, can 
be solved numerically to yield all g„^ s and V n ^ in terms of g\,\ 



which can be solved by y m ~ y Q m , with 9 = (2Q . We then 
have 



■p (\ m V n _}_rn—m 

1 n.m ~ " b — = ' 



(A3) and the flux (j) (n) scales with the layer index as 



B-2 



9"- 1 - 1 



m=0 



(A13) 



(A14) 
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Since 9 > 1, this is equivalent to 

4>(«)~(2Ce)" = i, (Ai5) 

so that the flux becomes asymptotically constant for n >• 1 . 
The function K (n), defined by 

2^—2 

K{n) = ^Go(x;x )Go(x;x ) = £ V„^ )S , (A16) 

xGn s=\ 

scales as (n) ~ and thus vanishes exponentially for large 
n. In the same limit, the function F («), related to K(n) and 
§ (n) through the equation 

J £F(t)K(n-t+l) = $(n), (A17) 
t=i 
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